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Abstract 



We explore a computational model of an incompressible fluid with 
^^.,^ a multi-phase field in three-dimensional Euclidean space. By inves- 

^ ■ tigating an incompressible fluid with a two-phase field geometrically, 

we reformulate the expression of the surface tension for the two-phase 
field found by Lafaurie, Nardone, Scardovelli, Zaleski and Zanetti (J. 
(S) ' Comp. Phys. 113 (1994) pp. 134-147) as a variational problem related 

r-«^ , to an infinite dimensional Lie group, the volume-preserving diffeomor- 

^D I phism. The variational principle to the action integral with the surface 

energy reproduces their Euler equation of the two-phase field with the 
surface tension. Since the surface energy of multiple interfaces even 
with singularities is not difficult to be evaluated in general and the 
p\ ' variational formulation works for every action integral, the new for- 

j^ . mulation enables us to extend their expression to that of a multi-phase 

(A'^-phase, N >2) flow and to obtain a novel Euler equation with the 
surface tension of the multi-phase field. The obtained Euler equation 
governs the equation of motion of the multi-phase field with different 
surface tension coefficients without any difficulties for the singular- 
ities at multiple junctions. In other words, we unify the theory of 
multi-phase fields which express low dimensional interface geometry 
and the theory of the incompressible fluid dynamics on the infinite di- 
mensional geometry as a variational problem. We apply the equation 
to the contact angle problems at triple junctions. We computed the 
fluid dynamics for a two-phase field with a wall numerically and show 



the numerical computational results that for given surface tension co- 
efficients, the contact angles are generated by the surface tension as 
results of balances of the kinematic energy and the surface energy. 
Keywords: multi-phase flowsurface tensionmultiple junction volume- 
preserving diffeomorphism 
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1 Introduction 

Recently, since the developments of both hardware and software in computer 
science enable us to simulate complex physical processes numerically, such 
computer simulations become more important from industrial viewpoints. 
Especially the computation of the incompressible multi-phase fluid dynamics 
has crucial roles in order to evaluate the behavior of several devices and 
materials in a micro- region, e.g., ink-jet printers, solved toners and so on. In 
the evaluation, it is strongly required that the fluid interfaces with multiple 
junctions are stably and naturally computed from these practical reasons. 

In this article, in order to handle the fluid interfaces with multiple junc- 
tions in a three dimensional micro-region, we investigate a surface tension 
of an incompressible multi-phase flow with multiple junctions as a numerical 
computational method under the assumption that the Reynolds number is 
not so large. In the investigation, we encounter many interesting mathemat- 
ical objects and results, which are associated with low dimensional interface 
geometry having singularities, and with the inflnite dimensional geometry of 
incompressible fluid dynamics. Further since even in a macroscopic theory, 
we introduce artiflcial intermediate regions in the material interfaces among 
different fluids or among a solid and fluids, the regions give a resolution of the 
singularities in the interfaces to provide extended Euler equations naturally. 
Thus even though we consider the multi-phase fluid model as a computational 
model, we believe that it must be connected with mathematical nature of real 
fluid phenomena as their description. We will mention the background, the 
motivation and the strategy of this study more precisely as follows. 

For a couple of decades, in order to represent the physical process with 
the interfaces of the multi-phase fluids, the computational schemes have been 
studied well. These schemes are mainly classifled into two types. The flrst 
type is based on the level-set method |42j discovered by H-K. Zhao, T. Chan, 
B. Merriman, S. Osher and L. Wang [l6l HS]. The second one is based on 



the phase-field theory, which was found by J. U. Brakbill, D. B. Kothe and 
C Zemach [TT], and B. Lafaurie, C Nardone, R. Scardovelli, S. Zaleski, and 
G. Zanetti [3^. The authors in Reference [31] called the scheme SURFER. 
Following them, there are many studies on the SURFER scheme, e.g., [71[T2| 
[211 references therein]. 

The level-set method is a computational method in which we describe a 
(hyper-) surface in terms of zeros of the level-set function, i.e., a real function 
whose value is a signed distance from the surface, such as q{x) in Section 
12.11 Using the scheme based upon the level-set method in the three dimen- 
sional Euclidean space, we can deal well with topology changes, geometrical 
objects with singularities, e.g., cusps, the multiple junctions of materials, 
and so on. However in the computation, we need to deal with the constraint 
conditions even for two-phase fluids |l6l HS]. A dynamical problem with 
constraint conditions is basically complicate and sometimes gives difficulties 
to find its solution since the constraint conditions sometimes generate an 
ill-posed problem in the optimization. In the numerical computation for in- 
compressible fluid, we must check the consistency between the incompressible 
condition and the constraint condition. The check generally requires a com- 
plicate implementation of the algorithm, and increases computational cost. 
Its failure sometimes makes the computation unstable, especially when we 
add some other physical conditions. Since instability disturbs the evaluation 
of a complex system as a model of a real device, it must be avoided. 

On the other hand, using the SURFER scheme [21], we can easily com- 
pute effects of the surface tension of a two-phase fluid in the Navier-Stokes 
equation. The phase fleld model is the model that we represent materials 
in terms of supports of smooth functions which roughly correspond to the 
partition of unity in pure mathematics [271 I p.272] as will be mentioned in 
Sections JH and O We call these functions "color functions" or "phase flelds" . 
The phase flelds have artiflcial intermediate regions which represent their 
interfaces approximately. In the SURFER scheme [31], the surface tension 
is given as a kind of stress force, or volume force due to the intermediate 
region. Hence the scheme makes the numerical computations of the surface 
tension stable. However it is not known how to consider a multi-phase [N- 
phase, N > 2) flow in their scheme. In Reference [11], the authors propose 
a method as an extension of the SURFER scheme [31] to the contact angle 
problem by imposing a constraint to flx its angle. In this article, we will 
generalize the SURFER scheme to multi-phase flow without any constraints. 

Nature must not impose any constraints even at such a triple junction. 



which is governed by a physical principle. If it is a Hamiltonian system, its 
determination must obey the minimal principle or the variational principle. 
We wish to find a theoretical framework in which we can consistently handle 
the incompressible fiows with interfaces including the surface tensions and 
the multiple junctions without any constraints. As the multiple junctions 
should be treated as singularities in a mathematical framework which are 
very difficult to be handled in general, it is hard to extend mathematical 
approaches for fiuid interface problems without a multiple junction [HI HO] to 
a theory for the problem with multiple junctions. Our purpose of this article 
is to find such a theoretical framework which enables us to solve the fiuid 
interface problems with multiple junctions numerically as an extension of the 
SURFER scheme. 

For the purpose, we employ the phase field model. The thickness of 
the actual intermediate region in the interface between a solid and a fiuid 
or between two fiuids is of atomic order and is basically negligible in the 
macroscopic theory. However the difference between zero and "the limit to 
zero" sometimes brings a crucial difference in physics and mathematics; for 
example, in the Sato hyperfunction theory, the delta function is regarded 
as a function in the boundary of the holomorphic functions [261 129], i.e., 

5{x) = lim — : I : — J = lim -. As mentioned above, the 
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phase field model has the artificial intermediate region which is controlled 
by a small parameter e and appears explicitly even as a macroscopic theory. 
We regard that it represents the effects coming from the actual intermediate 
region of materials. Namely, we regard that the stress force expression in 
the SURFER scheme is caused by the artificial intermediate region of the 
phase-fields and it represents well the surface effect coming from that of real 
materials. 

In order to extend the stress force expression of the two-phase flow to 
that of the multi-phase (iV-phase, N > 2) flow, we will flrst reformulate the 
SURFER scheme in the framework of the variational theory. In Reference 
[21] , a similar attempt was reported but unfortunately there were not precise 
derivations. Our investigations in Section Hj show that the surface tension 
expression of the SURFER scheme is derived as a momentum conservation 
in Noether's theorem [TOl El] and its derivation requires a generalization of 
the Laplace equation [30] as the Euler-Lagrange equation [3 HD], which is 
not trivial even for a static case. 

In order to deal with this problem in a dynamics case consistently, we 



should also consider the Euler equation in the framework of the variational 
principle. It is well-known that the incompressible fluid dynamics is geo- 
metrically interpreted as a variational problem of an infinite dimensional Lie 
group, related to diffeomorphism, due to V. I. Arnold [Ull], D. Ebin and J. 
Marsden [18], H. Omori [38] and so on. Following them, there are so many 
related works [SlEllSSlETlllIllSlllllin]. 

On the reformulation of the SURFER scheme [21] for the dynamical case, 
we introduce an action integral including the kinematic energy of the incom- 
pressible fluid and the surface energy. The variational method reproduces 
the governing equation in the SURFER scheme. 

After then, we extend the surface energy to that of multi-phase fields 
and add the energy term to the action integral. The variational principle of 
the action integral leads us to a novel expression of the surface tension and 
the extended Euler equation which we require. Using the extended Euler 
equation, we can deal with the surface tensions of the multi-phase flows, the 
multiple junctions of the of phase fields including singularities, the topology 
changes and so on. We can also compute a wall effect naturally and a contact 
angle problem. The computation of the governing equation is freed from any 
constraints, except the incompressible condition. 

In other words, in this article, we completely unify the theory of the 
multi-phase (iV-phase, N > 2) field and the theory of the incompressible fluid 
dynamics of Euler equation as an infinite dimensional geometrical problem. 

Contents are as follows: Section [2] is devoted to the preliminaries of the 
theory of surfaces in our Euclidean space from a low- dimensional differential 
geometrical viewpoint [M| [201 [T9] and Noether's theorem in the classical field 
theory [5], [TOl |23]. Section [3] reviews the derivation of the Euler equation to 
the incompressible fluid dynamics following the variational method for an 
infinite-dimensional Lie algebra based upon Reference [IB]. In Section HJ 
we reformulate the SURFER scheme [31j. There the Laplace equation for 
the surface tension and the Euler equation in Reference [31] are naturally 
obtained by the variational method in Propositions |8] and [TOl Section [5] is 
our main section in which we extend the theory in Reference [21] to that for 
a multi-phase flow and obtain the Euler equation with the surface tension 
of the multi-phase field in Theorem [2l The extended Euler equation for 
the multi-phase flow is derived from the variational principle of the action 
integral in Theorem [1] As a special case, we also derive the Euler equation 
to a two-phase field with wall effects in Theorem |3l In Section 6, using 
these methods in the computational fluid dynamics [13, [221 [21], we consider 



numerical computations of the contact angle problem of a two-phase field 
because the contact angle problem for the two-phase field circumscribed in 
a wall is the simplest non-trivial triple junction problem. By means of our 
scheme, for given surface tension coefficients, we show two examples of the 
numerical computations in which the contact angles automatically appeared 
without any geometrical constraints and any difficulties for the singularities 
at triple junctions. The computations were very stable. Precisely speaking, 
as far as we computed, the computations did not collapse for any boundary 
conditions and for any initial conditions. 

2 Mathematical Preliminaries 

2.1 Preliminary of surface theory 

In this subsection, we review the theory of surfaces from the viewpoint of 
low- dimensional differential geometry. The interface problems have been also 
studied for last three decades in pure mathematics, which are considered as a 
revision of the classical differential geometry [16] from a modern point of view 
[T71 [T9t [20| |3ll ST], e.g., generalizations of the Weierstrass-Ennpper theory 
of the minimal surfaces, isothermal surfaces, constant curvature surfaces, 
constant mean curvature surfaces, Willmore surfaces and so on. They are 
also closely connected with the harmonic map theory and the theory of the 
variational principle [191 120] ■ 

We consider a smooth surface S embedded in three dimensional Euclidean 
space ¥?. Let x = {x^,x^,x^) be of the Cartesian coordinate system and 
represent a point in E^, and let the surface 5* be locally expressed by a local 
parameter (s^, s^). We assume that the surface 5* is expressed by zeros of a 
real valued smooth function q over E'^, i.e., 

q{x) = 0, 

such that in the region whose |g| is sufficiently small (|g| < Et for a positive 
number Et > 0), \dq\ agrees with the infinitesimal length in the Euclidean 
space. Then dq means the normal co- vector field (one- form), i.e., for the 
tangent vector field Cq, := da '■= d/ds°' (a = 1, 2) of S, 

{da, dq) = over 5 = {x G E^ \q{x) = 0}. (2.1) 

Here (, ) means the pointwise pairing between the cotangent bundle and the 
tangent bundle of E^. The function q can be locally regarded as so-called the 
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level-set function [121 SS] . We could redefine the domain of q such that it is 
restricted to a tubular neighborhood T5 of S, 

Ts:={xeE^ I \q{x)\ < St}. 

Over Ts, q agrees with the level-set function of S. There we can naturally 
define a projection map n : Ts ^ S and then we can regard Ts as a fiber 
bundle over S, which is homeomorphic to the normal bundle A^^ — )■ S. How- 
ever the level-set function is defined as a signed distance function which is 
a global function over E^ as a continuous function |12] and thus it has no 
natural projective structure in general; for example, the level-set function L 
of a sphere with radius a is given by 

L{x\ x\ x^) = v/(xi)2 + (x2)2 + (x3)2 - a, 

which induces the natural projective (fiber) structure but the origin (0, 0, 0) 
in the sphere case. The level-set function has no projective structure at 
(0, 0, 0) in this case, and we can not define its differential there. In other 
words, the level-set function is not a global function over ¥? as a smooth 
function in general. 

When we use the strategy of the fiber bundle and its connection, we 
restrict ourselves to consider the function q in T5. Then the relation (12.11) 
and the parameter (si,S2) are naturally lifted to Ts as an inverse image of 

TT. 

Further for Cq := dq := d/dq, we have 

da{eq) = ^ r^^gC/? over 5". 

Here (F^ ) is the Weingarten map, which is a kind of a point- wise 2 x 2- matrix 
((^09)0^) [23 Chapter VII]. The eigenvalue of (Ff^) is the principal curvature, 
whereas a half of its trace tr(F^g)/2 is known as the mean curvature and its 
determinant det(F^g) means the Gauss curvature [271 Chapter VII]. 

Noting the relation, (e^, ds") = 6^ for a, /3 = 1, 2, the twice of the mean 
curvature, k, is given by, 

2_]9a{eq)ds°' = K over S. 

a 

Further noting the relation dqCqdq = 0, we obtain 

y da{eq)ds°' + dq{eq)dq = k over S. 



Due to the flatness of tlie Euclidean space, we identify e^ with Vg/| Vg| and 
then we have the following proposition. 

Proposition 1. The following relation holds at a point over S , 

div ( -r=-r I = K. 

Vivg|; 

For the case |Vg| = 1, using the Hodge star operator [5l |36] and the 
exterior derivative d, we also have an alternative expression *d*dq = k over 
the surface S. Here the Hodge star operator is * : IsP{Ts) — )• Is?~^{Ts) and 
the exterior derivative d : Ap(T5) -^ A^+^(rs) {du = Y.Li^i'^dx'), where 
Ap(Ts) is the set of smooth p-forms over Ts [5B] . 

Noting that as the left hand side of formula in Proposition [1] can be lifted 
to Ts, the formula plays an important role in References [HI [311 116] and in 
this article. 

2.2 Preliminary of Noether's theorem 

In this subsection, we review Noether's theorem in the variational method 
which appears in a computation of the energy-momentum tensor-field in the 
classical field theory [5l [TOl [23] . 

Let the set of i smooth real-valued functions over n-dimensional Eu- 
clidean space E" be denoted by C°°(E")®^, where n is mainly three. Let 
X = (a;^,x^, . . . ,x") be of the Cartesian coordinate system of E". We con- 
sider the functional / : C°°(E")®^ -^ M, 

/= / rf"x^(0„(x),9,0,(x)), (2.2) 

where J^ is a local functional, J^ : C°°(E'^)^% -^ A"(E")|^, 

J^ ■ (0a)a=i,.../U ^J^i(pa{x), <9j0a(a;))rf"x = J'(0„(x), <9i0a(a;), . . . , dn(l)a{x))d"-x 
= J'(0i(x), . . . , 4>£{x), di4>i{x), ..., dn(pe{x))d'^x 

and di := d/dx"^, {i = 1, ■ ■ ■ ,n). Then we obviously have the the following 
proposition. 



Proposition 2. For the functional I in /i2.^) over C°°(E")®^, the Euler- 
Lagrange equation coming from the variation with respect to (pa of{(f)b)b=i,...,i ^ 



C°°(E'')^\ i.e., 



51 



54>a{x) 



0, is given by 



SJ" 



E^^ 



67" 



Ha{x) j^ Sdi(f)a{x) 



0. 



(2.3) 



Using the equation (12.31) . we consider an effect of a small translation x 
to x' = X + Sx on the functional I. The following proposition is known as 
Noether's theorem which plays crucial roles in this article. 



Proposition 3. The functional derivative I with respect to 6xi is given by 

d^ [T] . (2.4) 






5x 



i=i 



X] T^n-^i'Pa 



a=l ^^^^^ 



If I is invariant for the translation, l[2.4^ gives the conservation of the mo- 
mentum. 

Proof. For the variation x' = x + 5x, the scalar function becomes 



<Pa{x') = Mx) + J2diMx)Sx' + 0(5. 

i=l 



X 



From the relations on the Jacobian and each component, 



dx' 



nrr •^— ^ 



dx 



we have 



i=l 



dx^ 
dx'"" 



5^^-d,5x'' + 0{5x^), 



[x 



{^) + ELi dj(l)a{x)5x^ dx 



dx' 



dx^ 



dx 



- + OiSx' 

I'l ^ 



di^a + E(9,9,-0„)(5X^' + 0(5 



X 



Then up to (5x^, we obtain 



(rx'F{(t)a{x),d[<Pa{x'))- / (r'xF{(t)a{x),d,<Pa{x)) 

E E T^^^0a(a;)'^a;* + E E T^— 5.5.0a (a;)'5a;^ + Yl ^^^^^'] ^"^ 

j=l a=l '^"- i,j=l a=l ^^°- j=l 

n i 



E" v»=l 



T.T.ji^^^^^-^ 



j = l a=\ 



Here we use the Euler-Lagrange equation (12. 3p and then we have (12. 4p . If we 
assume that / is invariant for the variation, it vanishes. D D 

3 Variational principle for incompressible fluid 
dynamics 

As we will derive the governing equation as the variational equation of an 
incompressible multi-phase flow with interfaces using the variational method, 
let us review the variational theory of the incompressible fluid to obtain the 
Euler equation following References [H H [TBI |25l EH |32l [37] . 

Let n be a smooth domain in '¥? . The incompressible fluid dynamics 
can be interpreted as a geometrical problem associated with an infinite di- 
mensional Lie group [H [HI [38] . It is related to the volume-preserving dif- 
feomorphism group SDiff(r2) as a subgroup of the diffeomorphism group 
Diff(i7). The diffeomorphism group Diff(i7) is generated by a smooth co- 
ordinate transformation of f2. The Lie algebras 5()iff(r2) = TeSDiff(f2) of 
SDiff (fi) and Oiff(f2) = TgDiff (f2) of Diff (il) are the infinite dimensional real 
vector spaces. The s^iff(f2) is a linear subspace of c)iff(f2). 

Following Ebin and Marsden [18], we consider the geometrical meaning 
of the action integral of an incompressible fluid, 

dt I d^x (-p\u\A . (3.1) 

Here T := (0,To) is a subset of the set of real numbers M, (x,t) is the 
Cartesian coordinate of the space-time Q x T, p is the density of the fluid 
which is constant in this section, and u = {u^,u'^,u^) is the velocity field of 
the fluid. 

10 



Geometrically speaking, a flow obeying the incompressible fluid dynamics 
is considered as a section of a principal bundle IFluid(r2 x T) over the absolute 
time axis T C M as its base space, 

SDifT(fi) > IFluid(r] x T) 

(3.2) 

T. 

The projection w is induced from the trivial fiber structure wq, : QxT ^ T, 
{{x,t) -^ t). In the classical (non-relativistic) mechanics, every point of 
space-time has a unique absolute time t G M, which is contrast to one in the 
relativistic theory. 

Due to the Weierstrass polynomial approximation theorem |1S], we can 
locally approximate a smooth function by a regular function. Let the set of 
smooth functions over Q be denoted by C°°{Q) and the set of the regular real 
functions by C^IQ) whose element can be expressed by the Taylor expansion 
in terms of local coordinates. 

The action of Diff (fi) on C^(fi) C C°°(l]) is given by 

e'-'^^f{x) = fix + su), 

for an element / G C^IQ), and small s > 0, where di := d/dx^ and we use the 
Einstein convention; when an index i appears twice, we sum over the index 
i. Thus the action e**"'^* is regarded as an element of Diff(n). 

As a frame bundle of the principal bundle IFluid(fi x T), we consider a 
vector bundle Coor(i7 x T) with infinite rank, 

C°°(fi) > Coor(fi X T) 



T. 

Since C°°(0) is regarded as a non-countably infinite dimensional linear space 
over R, we should regard Diff ($7) and SDiff (^2) as subgroups of an infinite 
dimensional general linear group if defined. 

More rigorously, we should consider the ILH space (inverse limit of Hilbert 
space) (or ILB space (inverse limit of Banach space)) introduced in Reference 
[38] by adding a certain topology to (a subspace of) C°°{Q x T), and then we 
also should regard Diff and SDifT as an ILH Lie group. However our purpose 
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is to obtain an extended Euler equation from a more practical viewpoint. 
Thus we formulate the theory primitively even though we give up to consider 
a general solution for a general initial condition. 

We consider smooth sections of Coor(r2 x T) and IFluid(f2 x T). Smooth 
sections of Coor(n x T) can be realized as C°°{Q x T). In the meaning of the 
Weierstrass polynomial approximation theorem |1S], an appropriate topology 
in C°°(r2 X T) makes C^{Q, x T) dense in C^{Q, x T) by restricting the region 
Q X T appropriately. Under the assumption, we also deal with a smooth 
section of IFluid(fi x T). 

Let us consider a coordinate function (7*(a;, t))i=i,2,3 G C'^{Q x T) such 
that 

—-f\x,t)=u\x,t), Y{x,t) = x' at teT, 

(JjL 

which means 

Y{x, t + 6t)=x' + u\x, t)6t + 0{5t^), 

for a small St. Here the addition is given as a Euclidean move in K^. As an 
inverse function of 7 = 'j{u, t), we could regard m as a function of 7 and t, 

u{x,t) = u{'y{x,t),t). 

Further we introduce a small quantity modeled on 5t ■ u^, 

f{x,t):=Y{x,t)-x\ (3.3) 

Then a section g of IFluid(i7 x T) at t G T can written by, 

g{t) = e^"^' G IFluid(fi xT) ^ SDiff(fi) C Diff(fi). (3.4) 

t 

Here we consider g as an element of SDiff (!]) and thus it satisfies the condition 
of the volume preserving, which appears as the constraint that the Jacobian, 

must preserve 1, i.e., the well-known condition that tr((9jU*) = div(u) must 

vanish, or —7— =0. 
at ox 
Following Reference [18], we reformulate the action integral (13. ip as "the 

energy functional" in the frame work of the harmonic map theory. In the 
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harmonic map theory [20] by considering a smooth map /i : M — t- G for a 
n-smooth base manifold M and its target group manifold G, "the energy 
functional" is given by 

E = l f ti {{h-^dh) * {h-^dh)) . (3.5) 

2 Jm 

Here * means the Hodge star operator, which is for * : TG ® A^(M) — )■ 
TG ® A"'^^(M) where AF{M) is the set of the smooth p-forms over M |36] . 
and TG® AP(M) is the set of the tangent bundle TG valued smooth p-forms 
over M [36j. The term "energy functional" in the harmonic map theory 
means that it is an invariance of the system and thus it sometimes differs 
from an actual energy in physics. 

Since in (13.21) . the base space T is one-dimensional and the target space 
IFluid(f2 X T)|i at t G T is the infinite dimensional space, "the energy func- 
tional" (13.51) in the harmonic map theory corresponds to the action integral 
i^freeb] which is defined by 



^-w -ILIJ^i^-'^^'^ ((-'"'*l»'"') (-'"4-"")) 



Here dx^{dj) := {dj^dx^) = 6^j is the natural pairing between TQ and T*Q. 
The trace in (13. 5 p corresponds to the integral over f2 with -^pd^x ■ dx^ dx^. 

The Hodge * operator acts on the element such as * (e~^ ^^dt-j^e' *) = 

( e~^ ^^ ^e^ ^'^ ) as the natural map from Oiff (i7) valued 1-form to 0-form. 
Further we assume that p is a constant function in this section. Then the 
action integral iSfree[7] obviously agrees with (13. ip . 

We investigate the functional derivative and the variational principle of 
this 5free[7]- Let us consider the variation, 

7^(0;, t') = 7^'(x,t') + 5^^{x,t'), and 7^(x,t') = i^{x,t') + 5^^{x,t'), 

where we implicitly assume that 5''^^ is proportional to the Dirac 5 function, 
5{t' — t), for some t and Sj^ vanishes at dfl. As we have concerns only 
for local effects or differential equations, we implicitly assume that we can 
neglect the boundary effect arising from dQ on the variational equation. If 
one needs the boundary effect, he would follow the study of ShkoUer |13] . 
Further one could use the language of the sheaf theory to describe the local 
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effects [25]. As we are concerned only witli differential equation and thus 
our theory is completely local except Section El we could deal with germs 
of related bundles [2] as in Reference [31] , which is also naturally connected 
with a computational method of fluid dynamics [35] . 

Let us consider the extremal point of the action integral (13.11) following the 
variational principle. Noting that d^jdx = 1, the above Jacobian becomes 



dx dx 



:i + dk5^') + Oii5^y). 



Since we employ the projection method, we firstly consider a variation in 
^iff(fi) rather than 50iff(fi). For the variation, the action integral iSfree[7] 
with (13.41) becomes 



-j^dtj^^d'x ■ dx^ ® dx^ (^57^1 {P9''jt9) + ^l'd,lp\u\'^ 



Now we have the following proposition. 

Proposition 4. Using the above definitions, the variational principle in 
SDiff(l]), 



= 0, 

SDiff(0)|t 



6-f{x,t) 
is reduced to the Euler equation, 

—pu' + u^djpu' + d,p = 0, (3.6) 

where p comes from the projection from TDiff(r2)|sDiff(Q) — ^ TSDiff(i7). 

Proof. Basically we leave the rigorous proof and especially the derivation of p 
to [H HHj. The existence of p was investigated well in Appendix of Reference 
[T8] as the Hodge decomposition [5], [36|. (See also the following Remark 
[H) Except the derivation of p, we use the above relations and the following 
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relations, 



(9 . .\ 

—pu' + u^djpu' 1 di 

Then we obtain the Euler equation. D D 

Remark 1. The Euler equation was obtained by the simple variational prin- 
ciple. Physically speaking, the conservation of the momentum in the sense 
of Noether's theorem [TOl |23] led to the Euler equation. However, we could 
introduce the pressure pl term as the Lagrange multiplier of the constraint 
of the volume preserving. In the case, instead of iSfree, we deal with 



'^free.p = '^free + dt pL{x,t)—d^X. 



Then noting the term coming from the Jacobian, the relation. 



= 0, 

SDiff(f7)|t 



6-f{x,t) 

is reduced to the Euler equation, 

d 1 

—pu' + uWjpu' + di{pL + -p\u\^) = 0. 

As the pressure is determined by the (divergence free) condition of u, we 
renormalize [281 (25)], 

1 , ,2 
P-=PL + -p\u\ . 

More rigorous arguments are left to References fi8[ [38j and physically inter- 
pretations are, e.g., in References [HI [HI [251 133 SB 1^ - 
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We give a comment on the projection from TDiff(r2)|sDiflf(f7) -^ TSDiff(O) 
in fl3.6p . which is known as the projection method. First we note that the 
divergence free condition div(M) = simphfies the Euler equation fl3.6p . 

Dti f)ii^ 1 

p— + Vp = 0, — + u^d,u' + -d,v = 0. 
Dt at p 

As mention in Section |6], in the difference equation we have a natural in- 
terpretation of the projection method [I3]. We, thus, regard Du/Dt in 
TDiff(n)|sDiff(n) as hm^i^o ^^^^±^p^ for u{t + 5t) := u{t + 5t,^{t + 5t)) e 
()iff(f2) and u{t) := u{t,^{t)) E sDiff(fi), i.e., div(ti(t)) = by considering 
TDiff(fi) at the unit e of Diff ((]) up to fe^, as we did in (Q and (jSl]). In 
order to find the deformation ^"(t + 5t) in sOiff(r2) by a natural projection 
from ()iff(r2) to sDiff(r2) [14, ,p.36], we decompose u{t + 5t) into u^t + 5t) and 
u^{t + 6t) such that diU^%t + 5t) := diu\t + 5t). Then u\\t + 5t) := u{t + 5t)- 
u^(t + 6t) belongs to S()iff(fi). Thus the pressure p is determined by [1^ 

diu\t + 6t) + 6td^-d^p = 0. (3.7) 

P 



In other words, since ^"(t + 6t) = u^(t + 5t) + 5t-dip belongs to sDiff(r2), the 

deformation of M"*(t + (5t) —u^{t) which gives Du^Dt and the Euler equation 
(13.61) is the deformation in IFluid(fi x T). After taking the continuous limit 
6t — > 0, the equation for the pressure (13.71) can be written as [13], 

idiu')idju')+di-dip = 0, 

by noting the relations [dt,di] = and div(M(t)) = 0, i.e., diU^it + 5t) = 
di[u'{t) +§-^u\t)5t +u^{t)djv}{t)5t\+0{6t'^). The Poisson equation with dSSJ 
guarantees the divergence free condition. Hence the pressure p in the incom- 
pressible fluid is determined geometrically. 



4 Reformulation of Surface tension as a min- 
imal surface energy 

In this section we reformulate the SURFER scheme [31] following the varia- 
tional principle and the arguments of previous sections. 
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4.1 Analytic expression of surface area 

We first should note that in general, the higher dimensional generalized func- 
tion like the Dirac delta function has some difficulties in its definition |18]. 
For the difficulties, in the Sato hyperfunctions theory [26] , the sheaf theory 
and the cohomology theory are necessary to the descriptions of the higher 
dimensional generalized functions, which are too abstract to be applied to 
a problem with an arbitrary geometrical setting. Even for the generalized 
function in the framework of Schwartz distribution theory, we should pay 
attentions on its treatment. However since the surface S in this article is a 
hypersurface and its codimension is one, the situation makes the problems 
much easier. 

We assume that the smooth surface S is orientable and compact such that 
we could define its inner side and outer side. In other words, there is a three 
dimensional subspace (a manifold with boundary) B such that its boundary 
dB agrees with S and B is equal to the inner side of S with S itself. Then 
we consider a generalized function 6 over f2 C E^ such that it vanishes over 
the complement B'^ = Q\ B and is unity for the interior B° := B \ dB; 9 is 
known as a characteristic function of B. 

We consider the global function 9{x) and its derivative d9{x) in the sense 
of the generalized function, which is given by 

d9{x) = J2 di9{x)dx' = dq9{x)dq. 

i 

Here we use the notations in Section 12.11 Using the nabla symbol V^ = 
{di9{x))i=i^2,3, \V9\d^x is interpreted as 

|V^|rf^x= \{*d9)Adq\. 

Here due to the Hodge star operation * : AP(r2) — )■ A'^^P(f2), *d9 = edg9ds^ A 
ds'^ where e is the Jacobian between the coordinate systems {ds^, ds"^, dq) and 
((ix^, (ix^, dx^). Then we have the following proposition; 

Proposition 5. // the integral, 

A:= f \V9\d^x= f \{*d9)Adq\, 
Jn Jn 

is finite, A agrees with the area of the surface S . 

It should be noted that due to the codimension oi S C Q, we have used 
the fact that the Dirac 6 function along g G T^ is the integrable function 
whose integral is the Heaviside function. This fact is a key of this approach. 
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4.2 Quasi-characteristic function for surface area 

For the later convenience, we introduce a support of a function over Q, which 
is denoted by "supp", i.e., for a function g over Q, its support is defined by 

supp(5f) = {x en \ g{x) y^ 0}, 

where " " " means the closure as the topological space Q. 

One of our purposes is to express the surface S by means of numerical 
methods, approximately. Since it is difficult to deal with the generalized 
function 6* in a discrete system hke the structure lattice [15], we introduce a 
smooth function ^ over Q as a quasi-characteristic function which approxi- 
mates the function 9 [HI [31], 

( for X e B^f]{x E n I \q{x)\ < eg/2}^ 

^{x)=< 1 ioi X eBf]{xen\\q{x)\<e^/2Y 

[^ monotonically increasing in q(x) otherwise. 

(4.1) 
We note that along the line of dq for q G (—eg/2, eg/2), ^ is a monotonically 
increasing function which interpolates between and 1. We now implicitly 
assume that eg is much smaller than et defined in Section [271] so that support 
of |V^| is in the tubular neighborhood T5. However after formulating the 
theory, we extend the geometrical setting in Section 12.11 to more general ones 
which include singularities; there Et might lose its mathematical meaning but 
eg survives as a control parameter which governs the system. For example, 
as in Reference [31], we can also deal with a topology change well. 

By letting ^'^{x) := 1 — ^{x), ^'^ and ^ are regarded as the partition of 
unity [271 I P-272], or 

e(x)+r(x) = i. 

We call these ^ and ^^ "color functions" or "phase fields" in the following 
sections. We have an approximation of the area of the surface S by the 
following proposition. 

Proposition 6. Depending upon eg, we define the integral, 

^g := / \V^\d^x, 



Jo, 
and then the following inequality holds, 

l^g - ^1 < eg • ^. 
18 



Here we note that A^ is regarded as the approximation of the area Aoi S 
controlled by e^. In other words, we use e^ as the parameter which controls the 
difference between the characteristic function 6 and the quasi-characteristic 
function ^ in the phase field model [HI [31] . 

Let us consider its extremal point following the variational principle in a 
purely geometrical sense. 

Proposition 7. For sufficiently small e^, we have 



= k{x), 
where x ^ S or q = 0. 
Proof. Noting the facts that d^/dq < at g = and 

|vei = Vve ■ ve, 

Proposition [2] and the equality in Proposition [1] show the relation. D D 

In the vicinity of S, q in Section 12.11 could be identified with the level-set 
function and the authors in References ^G] J35| also used this relation. Since 
all of geometrical quantities on S are lifted to Ts as the inverse image of %, 
the relation in Proposition [7] is also defined over (supp(|V^|))° C Ts and we 
redefine the k, by the relation from here. 

4.3 Statics 

Let us consider physical problems as we finish the geometrical setting. Before 
we consider dynamics of the phase field flow, we consider a statical surface 
problem. Let a be the surface tension coefficient between two fluids corre- 
sponding to ^ and ^'^. Now let us call ^ and ^^ "color functions" or "phase 
fields" . More precisely, we say that a color function with individual physical 
parameters is a phase field. The surface energy S := a A is, then, approxi- 
mately given by 

£two ■.= aA^ = (T f \V^\d^x. (4.2) 

Jn 
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As a statical mechanical problem, we consider the variational method of this 
system following Section 12.21 

Since a statical surface phenomenon is caused by the difference of the 
pressure of each material, we now consider a free energy functional |33] , 

^two := / (a|V^| - {pi^ + P2n) d^x, (4.3) 

Jn 

where Pa (a = 1, 2) is the proper pressure of each material. 

Proposition 8. The variational problem with respect to C,, SJ-'two/^C, = 0, 
reproduces the Laplace equation [30t Chap. 7], 

{pi - P2) - (Jk{x) = 0, X G (suppd Ve|))°. (4.4) 

Proof. As in Proposition [21 direct computations give the relation. D D 

This proposition implies that the functional J-'two is natural. The solu- 
tions of (14. 4p are given by the constant mean curvature surfaces studied in 
References [HI |20l ST] . 

Furthermore we also have another static equation, whose relation to the 
Laplace equation (14. 4p is written in Remark |H 

Proposition 9. For every point x eVL, the variation principle, 5J^x.vfo/5x^ = 
0, gives 

" I E ^'nff - E 3,^ ] - (Pi - P^)9.f = 0. (4.5) 



or 



where 



djTij{x) - {pi - P2)dii{x) = 0, (4.6) 



^(^):=^(/-^®^)|V^|(x). 

Proof. We are, now, concerned with the variation x ^ x + 6x for every point 
X G fi. We apply Proposition [3] to this case, i.e., 



5x^ 



—a 



^•i^«i-^'l^-«w-s4^i^«i 



[x) + (pi -p'2)dii{x) 



by using (14. 4p as its Euler-Lagrange equation (12.30 . Further for x (supp(| V^|))'' 
its Euler-Lagrange equation ( 12.30 gives a trivial relation, i.e., "0 = 0". Then 
we have (jMl). □ □ 
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Remark 2. It is worthwhile noting that fl4.5p and fl4.6p are defined over Q 
rather than (supp(|V^|))° because due to the relation, 

m\ < mi 

even at the point at which denominators in the first term in (14. 5p vanish, the 
first term is well-defined and vanishes. 

Hence fl4.5p and (14.60 could be regarded as an extension of the defined 
region of (14.41) to Q and thus (14.51) and (14. 6 p have the advantage over (14. 4p . 
The extension makes the handling of the surface tension much easier. 

Remark 3. In the statical mechanics, there appears a force diTij, which 
agrees with one in (33) and (34) in Reference [31] and (2.11) in Reference 
[21]. We should note that in Reference [23], it was also stated that this 
term is derived from the momentum conservation however there was not its 
derivation in detail. The derivation of the above r needs the Euler-Lagrange 
equation (12. 3p . which corresponds to the Laplace equation (14. 4p in this case, 
when we apply Proposition [3] to this system, though these objects did not 
appear in Reference j24j . 



Remark 4. In this remark, we comment on the identity between (14. 4p and 
(14.60 . Comparing these, we have the identity, 

diTij = adji ■ K, 

which is, of course, obtained from the primitive computations. It implies 
that (14. 6 P can be derived from the Laplace equation (14.40 with this relation. 
However it is worthwhile noting that both come from the variational prin- 
ciple in this article. In fact, when we handle multiple junctions, we need 
a generalization of the Laplace equations over there like (15.71) . which is not 
easily obtained by taking the primitive approach. Further the derivations 
from the variational principle show their geometrical meaning in the sense of 
References [HEIIIO]. 



4.4 Dynamics 

Now we investigate the dynamics of the two-phase field. There are two dif- 
ferent liquids which are expressed by phase fields ^ and ^^ respectively. We 
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assume that they obey the incompressible fluid dynamics. As in the previ- 
ous section, we consider the action of the volume-preserving diffeomorphism 
group SDiff (n) on the color functions ^ and i^'^. We extend the domain of 
^ and ^'^ to n X T and they are smooth sections of Coor(r2 x T). For the 
given t, we will regard ^ and C,'^ as functions of 7* in the previous section, 
i.e., C, = C,{j{x,t)). For example, the density of the fluid is expressed by the 
relation, 

P = PiC + P2i 

for constant proper densities pi and p2 of the individual liquids. The density 
p, now, differs from a constant function over fi x T in general. 

We consider the action integral 5two including the surface energy, 

5two[7] = jj^j^ (^/'l^l' - ^l^^l + (M +P2r)) d^x. (4.7) 

The ratio between p and a determines the ratio between the contributions of 
the kinematic part and the potential (or surface energy) part in the dynamics 
of the fluid. Since the integrand in (14.71) contains no dS,/dt term, we obtain 
the same terms in the variational calculations from the second and the third 
term in (14. 7p as those in (14. 4p and (14.61) in the static case even if we regard 
?T, as 4 and x^ as t in Section 12.21 By applying Proposition [2] to this system, 
we have the following proposition as the Euler-Lagrange equation for ^. 

Lemma 1. The function derivative of Stwo with respect to ^ gives 

-ipi-p2)\u{x,t)\^ + {pi-p2)~at^{x,t) = 0, X e (supp(|Ve|))°, (4.8) 

up to the volume preserving condition. 

This could be interpreted as a generalization of the Laplace equation (14.41) 
as in the following remark. 

Remark 5. Here we give some comments on the generalized Laplace equa- 
tion (14. 8 p up to the volume preserving condition. This relation (14.80 does 
not look invariant for Galileo's transformation, m — )■ m + mq for a constant 
velocity uq. However for the simplest problem of Galileo's boost, i.e., static 
state on a system with a constant velocity Uq, the equation (14.81) gives 

-ipi-p2)\uo\^ + {pi-p2)-(TK{x,t) = 0, X G (supp(|Ve|))°, (4.9) 
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which might differ from the Laplace equation fl4.4p . However for the boost, 
we should transform Pa into 

Pa ■= Pa + -Pa\Uo\^ . (4.10) 



Then the above equation of pa agrees with the static one (14.41) . In other 
words (I4.10p makes our theory invariant for the Gaililio's transformation. 

For a more general case, we should regard pa as a function over Q x T 
rather than a constant number due to the volume preserving condition. These 
values are contained in the pressure as mentioned in (14.121) . The statement 
"up to the volume preserving condition" has the meaning in this sense. In 
fact, in the numerical computation, these individual pressures paS are not so 
important as we see in Remark [61 Due to the constraint of the incompressible 
(volume-preserving) condition, the pressure p is determined as mentioned in 
Remark [H There are no contradictions with the Galileo's transformation and 
SDiff(r2)-action. 

We consider the infinitesimal action of SDiff (fi) around its identity. As 
did in Section [31 we apply the variational method to this system in order to 
obtain the Euler equation with the surface tension. 

Proposition 10. For every {x,t) G QxT, the variational principle, SStwo/^Yi^yt) 
0, gives the equation of motion, or the Euler equation with the surface ten- 
sion, 

-^ + -[229^^^'229^^^)+9^P = ^- (4-11) 

Here p is also the pressure coming from the effect of the volume-preserving. 

Proof. The measure d^x is regarded as 7— rf x with 7— = 1. Noting --t— = 0, 

ox ox dt ox 

the proof in Proposition [H and Remark [H provide the kinematic part with 

pressure term and Proposition [91 gives the remainder. In this proof, the total 

pressure p is defined in Remark [6l D D 

Remark 6. More rigorous speaking, as we did in Remark [H we also renor- 
malize the pressure, 

P = Pl + ttPI^P + Pi^ + P2C 



2' 

1 1 

PL + -(Pl - P2)^\u\^ + (pi - P2)^ + ^P2\u\^ + P2- 
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(4.12) 



As in Section 12. 2[ the third term in fl4.1ip includes the effects from p^'s via 
the generahzed Laplace equation fl4.8p as the Euler-Lagrange equation ( I2.3p . 

Remark 7. 1. The equation of motion ( 14. lip is the same as (24) in Ref- 
erence pi] basically. We emphasize that it is reproduced by the varia- 
tional principle. 

2. As in Reference [21], in our framework, we can deal with the topology 
changes and the singularities which are controlled by the parameter e^. 
The above dynamics is well-defined as a field equation provided that 
e^ is finite. If needs, one can evaluate its extrapolation for vanishing of 

3. In general, eg is not constant for the time development. Due to the 
equation of motion, it changes. At least, in numerical computation, 
the numerical diffusion makes the intermediate region wider in gen- 
eral. However even when the time passes but we regard it as a small 
parameter, the approximation is justified. 

4. Since from Remark [21 the surface tension is defined over Q, the Euler 
equation is defined over Q without any assumptions. 

5. It should be noted that the surface force is not difficult to be computed 
as in Reference [21] but there sometimes appear so-called parasite cur- 
rent problems in the computations even though we will not touch the 
problem in this article. 

5 Multi-phase flow with multiple junctions 

In this section, we extend the SURFER scheme [31] of two-phase flow to 
multi-phase (iV-phase, N >2) flow. 

5.1 Geometry of color functions 

In order to extend the geometry of the color functions in the previous section, 
we introduce several geometrical tools. First let us deflne a geometrical object 
similar to smooth d-manifold with boundary. Here we note that (i-manifold 
means (i-dimensional manifold, and d-manifold with boundary means that its 
interior is a d- manifold and its boundary is a (rf — l)-dimensional manifold. 
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We distinguish a smooth (differential) manifold from a topological manifold 
here. 

When we consider mult i-j unctions in E^, we encounter a geometrical ob- 
ject with smooth "boundaries" whose dimensions are two, one and zero even 
though it is regarded as a topological 3-manifold with boundary. 

Definition 1. We say that a path- connected topological d-manifold with 
boundary V is a path- connected interior smooth d-manifold if V satisfies 
the followings: 

1. The interior V° is a path- connected smooth d-manifold, and 

2. V has finite path- connected suhspaces Va, (a = 1, ■ ■ ■ ,£) such that 

(a) V \V° is decomposed by Va, i.e., 



v\v° = Wv, 



a; 



(b) Each Va is a path- connected smooth k-manifold in Q {k < d). 

We say that Va is a singular-boundary ofV and let their union V\V° denoted 
byd,^^V:=V\V°. 

Here the disjoint union is denoted by ]J, i.e., for subsets A and B of Vt, 
AY[B := A{jB ii Ap\B = {h. 

By letting F^") := V and V^^"^ := {Va C V \ dim \4 < k], and by picking 
up an appropriate path-connected part V^^' C V^''^ each A;, we can find a 
natural stratified structure, 

y{n) ^ y{n-l) 3 ... 3 y{2) 3 y{l) 3 y{0) ^ 

which is known as a stratified submanifold in the singularity theory [2]. 

In terms of path-connected interior smooth (i-manifolds, we express sub- 
regions corresponding to materials in a regions i7 C E^ as extensions of B 
and -B^ in Section I^?T1 

Definition 2. For a smooth domain Vt C E^, we say that N path- connected 
interior smooth 3-manifolds {Ba}a=o.--,N-i (ii"^ colored decomposition ofQ if 
{Ba}a=o,-,N-i satisfy the followings. ■ 
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1. every Ba is a closed subset in Vt, 
2- ^ = [Ja=Q,-,N-iBa, and 

3. fi\(a<,fiani?,) = U.=o,...,iv-i^:- 

Roughly speaking, each Ba corresponds to a material in f2; Definition |2]1. 
means that Ba is surrounded by singular boundary or the boundary of fi, 2. 
implies that there is no "vacuum" in Vt and 3. guarantees that the interiors 
of these materials don't overlap. 

In general, for a 7^ 6, 5^ fl i?;, is a singular geometrical object if it is 
not the empty set. Singularity basically makes its treatment difficult in 
mathematics. In order to avoid such difficulties, we introduce color functions 
ia{.x) (a = 0, 1, 2, ■ ■ ■ , A^ — 1) over a region f2, which are modeled on ^ and ^^ 
as in Section 14. ![ are controlled by a small parameter e^ > and approximate 
the characteristic functions over Ba- 

To define color functions ia{,x) (a = 0, 1, 2, ■ ■ ■ , A^ — 1), we introduce 
another geometrical object, e-tuhular neighborhood in ¥?: 

Definition 3. For a closed suhspace U G Q and a positive number e, e- 
tubular neighborhood Tu,e of U is defined by 

Tu,e:={xen I dist(x,f/) < |}, 

where dist(x, U) is the distance between x and U inK^. 

We assume that each Tg^^^_^^Ba,e has a fiber structure over dsmgBa as topo- 
logical manifolds as mentioned in Section 12.11 Using the e-tubular neighbor- 
hood, we define e^-controUed color functions. 

Definition 4. We say that N smooth non-negative functions {C,a}a=o,--,N-i 
over fi C E^ are t£_- controlled color functions associated with a colored de- 
composition {Ba}a=i),--,N-i of Q, if they satisfy the followings: 

1. ^a belongs to C°°{VL) and for x eVL, 

a=0,l---,Af-l 
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2. For every Ma := supp(^a) and La := supp(l— ^q), (a = 0, 1, ■ ■ ■ , A^— 1), 

(a) Ba^Ma, 

(b) K g Ba, 

(C) iMa\LirCT9,^^B.,,^, 
(d) {Ma\Liy^d,^^Ba. 

3. For X G {Ma \ La), we define the smooth function qa by 

(^\ ^ / dist(x, dsingBa), X e {Ma \ Ba), 

|^—dist(x, casing -Ba), othtrwise. 

Then for the flow exp(— t^) on C°°(i7)|(7v/„\Lg), ^a monotonically in- 
creases along tEUcM.atxE {Ma \ L^). 

When {Ma \ L^f = Td^,^^Ba,e^ for every a, {U}a=o,-,N-i are called proper e^- 
controlled color functions associated with the colored decomposition of Q G 
K^, {Ba}a=o,---,N~-i or merely proper. 

The functions .^a's are, geometrically, the partition of unity [27i I p. 272] 
and a quasi-characteristic function, roughly speaking, which is equal to 1 in 
the far inner side of Ba, vanishes at the far outer side of Ba and monotonically 
behaves in the artificial intermediate region. Noting that the flow exp(— t^) 
corresponds to the flow from the outer side to the inner side, ^a decreases 
from the inner side to the outer side. 

From here, let us go on to use the notations Ba, Ma, La, and ^a in Defi- 
nition m Further for later convenience, we employ the following assumptions 
which are not essential in our theory but make the arguments simpler. 

Assumptions 1. We assume the following: 

1. The colored decomposition {Ba}a=o,--,N-i of Q and e^ satisfy the con- 
dition that every La is not the empty set. 

This assumption means that the singularities that we consider can be 
resolved by the above procedure. Since e^ can be small enough, this 
assumption does not have crucial effects on our theory. 
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2. The colored decomposition {Ba}a=o,--,N-i of Q and eg satisfy the rela- 
tion, 

\af^b;a,bf^O 

and every intersection Ba H Bq perpendicularly intersects with dVt. 

This describes the asymptotic behavior of the materials. For example 
Mq will be assigned to a wall in Section O This assumption is nei- 
ther so essential in this model but makes the arguments easy of the 
boundary effect. As mentioned in Section [3l we neglect the boundary 
effect because we are concerned only with a local theory or differential 
equations. If one wishes to remove this assumption, he could consider 
smaller region Vt' dVt after formulates the problems in VL. 

3. The volume of every Ba, the area of every dsmgBa, and the length defined 
over every one- dimensional object in dgi^gBa are finite. 

As our theory is basically local, this assumption is not essential, either. 

Under the assumptions, we fix colored decomposition {Ba}a=Q,---,N-i and 
e^-controUed color functions {^a}a=o,--,N-i- 

As mentioned in the previous section, we have an approximate description 
of the area of dgmgBa- 

Proposition 11. By letting the area of dsingBa be Aa, the integral 

Jn 



approximates Aa by 



\Aa - Al < e^A- 



Here we notice that Mah := M^ f) ^fe (a, 6 = 0, 1, 2, ■ ■ ■ , A^ - 1, a ^ 6) 
means the intermediate region whose interior is a 3-manifold. Similarly we 
define Mabc := M^ H ^fe fl ^c (a, &, c = 0, 1, 2, ■ ■ ■ , A^ - 1; a 7^ 6, c; 6 7^ c) and 
so on. Since the relation, IJ Ma = ^, holds, we look on the intersections of 
Ma as an approximation of the intersections of Ba which is parameterized 
by eg. Even though there exist some singular geometrical objects in {Ba} 
[2], we can avoid its difficulties due to finiteness of eg. (We expect that the 
computational result of a physical process might have weak dependence on 
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e^ which is small enough. More precisely the actual value is obtained by 
the extrapolation of e^ = for series results of different e^'s approaching to 
e^ = 0.) 

5.2 Surface energy 

Let us define the surface energy inexact by 

a>b 

where aab is the surface tension coefficient {aab > 0, aab = crta) between the 
materials corresponding to Ba and Bb, and Area(?7) is the area of an interior 
smooth 2-manifold U. 

We have an approximation of the surface energy inexact by the following 
proposition. 

Proposition 12. The free energy, 

^^^^ = $^a,, / d'^x V\'^Ux)\\VUx)\{U^)+U^)), (5.1) 

a>b -^^ 

has a positive number M such that 

Proof. For a ^ b, Baf]Bi) consists of the union of some interior smooth 
2-manifolds. Their singular-boundary parts dsing{Baf]Bb) = {Va}a£iab ^^^ 
union of some smooth 1-manifolds and smooth 0- manifolds. Thus {Va}aeiab 
has no effect on the surface energy ^exait because they are measureless. 
Over the subspace, 

M^" := {x e Mab I Ux) + Ux) = 1}°, (5.2) 

and for a positive number £, we have identities, 

\^Ux)\{Ux) + Ux)Y = \^Ux)\{Ux) + Ux)Y 

= VWUx)\mMiUx) + Ux)Y. 

The sum of the integrals over M^™'' dominates £^^^ if e^ is sufficiently small. 
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We evaluate the remainder. For example, for different a, b and c, the part 
in S^^") coming from 

^rr := {^ e Mate I U^) + U^) + U^) = 1}° (5.4) 

is order of e^^. Namely we have 

d'x ^/WU^^WU^liU^) + U^)) - Length(5„ nB^nB, 
<e^^Length{Ban BbHB,), 



Mabc 

2t 



where Length(C) is the length of a curve C. Thus we find a number M 
satisfying the inequality. D D 

Remark 8. 1. M is bound by 

M < max(a„5) I ^ (Area(5„ n B^) + e^Length (9sing(5« n Bt))) + Kt^^ , , 

\a<h 

where K is the number of isolated points in all of singular-boundary 
parts of {Sq}. 

2. It should be noted that 8^^^ becomes the surface energy of the system 
exactly when e^ vanishes. 

3. Using the identities ( 15. 3p . we can also approximate £^*^^-' by 

^aab I d?x \VU^)\{Ux)+Ux)Y, 

a>b -^^ 

using a positive number £. In such a way, there are so many variants 
which, approximately, represent the surface energy in terms of ^^'s. 

5.3 Statics 

Let us consider the statics of the multi-phase fields. In the above arguments 
in this section, we have given the geometrical objects, first, and defined the 
functions ^a, functional energy £^^^ and so on. In this subsection on the 
static mechanics of the multi-fields, we consider the deformation of these 
geometrical objects and determine a configuration whose corresponds to an 
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extremal point of the functional, i.e., J^mui in the following proposition. In 
other words, we derive the Euler-Lagrange equation which governs the ex- 
tremal point of the functional and characterizes the configuration of Ma, La 
and approximately Ba for every a = 0, ■ ■ ■ ,N — 1. 
Let us introduce the proper pressure 



X) 



where pa is a certain pressure of each material. 



(5.5) 



Proposition 13. The Euler-Lagrange equation of the static free energy in- 
tegral, 



^. 



inul 



J2 ^abVl^^U^WUxMai^) + U^)) -Pp]dr 






,a>b 



with respect to ^a, i-S-., SJ^mxA/^ia = 0, is given as follows: 
L For a point x G M^^p^p of / fOI) . 

(Pa-Pb) -(^abf^aix) =0, 

where 



(5.6) 



Ka '■ — di 






2. For a point x E M^^^ of ^5^^, 



{Pa -Pb- Pc) - K,abc{x) = 0, 



(5.7) 



where 



kabc := (TbcVWU^^^WU^l - ^afe\/|Vea(x)||Vefe(a;)| - ^ac V^V^MM^ 



,^3 - Uab^/\V^\{^a + 6) + (^ac^/\VUi^a + ^c 



(5.8) 



Proof. For a point x G Mab^^°^ of (15. 2p . we have ^a{x) + C,b{x) = 1, and thus 
the Euler-Lagrange equation (12.31) leads (15.61) . 
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Similarly for a point x G Mq^c''™'' of (15.41) . we have ia{x)+^i,{x)+^c{.x) = 1, 
and thus the concerned terms of the integrand in the energy functional are 
given by 

■ ■ ■ + VWu^WU^li^a + 6) + VWu^Wu^Ku + a) .g g. 
+ v/|va^)iiva^)i(i-^„) + ---. 

The Euler-Lagrange equation (12. 3p gives (15. 7p . D D 

Remark 9. 1. It is noticed that (\5.Q\\ agrees with the Laplace equation 
(14.41) and thus we also reproduce the Laplace equation locally. 

2. (15. 7p could be regarded as another generalization of the Laplace equa- 
tion though M^l°^ does not contribute to the surface energy when e^ 
vanishes and has a negligible effect even for a finite e^ if eg is sufficiently 
small. Indeed, (15. 7p does not appear in the theory of surface tension 
[30] . However (15. 7p is necessary and plays a role to guarantee the sta- 
bility in the numerical computations and to preserve the consistency in 
numerical approach with finite intermediate regions for e^ 7^ 0. 

3. Similarly we have similar equations for a higher intersection regions. 

As a generalization of (14. 5 p we immediately have the following. 

Proposition 14. For every point x G Q, the variational principle, 6J^raui/Sx'^ = 
0, gives 

^^PP - J2 ^«.^ [^* {VWUl^li^a + 6)) 
a>b 

\VNU )■ 

Proof. It is the same as Proposition |9l which essentially comes from Propo- 
sition El D D 

Remark 10. In Proposition [TH we can apply the equation without any 
classification of geometry like (15. 2p and (15. 4p . It is also noted that (I5.10p is 
globally defined over Vt as mentioned in Remark |2J 
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5.4 Dynamics 

Using these equations, let us consider the dynamics of the multi-phase flow. 
We extend the colored-decomposition of Q and the e^-controlled color func- 
tions of {^a}a=o,---,N-i to those of ^2 X T and C°°{Q x T) using another flber 
structure of Coor(r2 x T). Mathematically speaking, since our space-time 
is a trivial bundle Q x T and has the flber structure Q x {ta,tb) — )■ Q for a 
small interval (ta,tb) due to the integrability, we can consider the pull-back 
of the map ^^ : i7 — ;■ M. If we consider a global behavior of C,a with respect 
to time t, we should pay more attentions on the Lagrange picture 7(x, t) and 
the integrability. However as our theory is local, we can regard (ta,tb) as T 
with an inflnitesimal interval. 

Thus ^a is redeflned as ^a '■= ^a{l{x,t)) for (x,t) G Q x T and it is 
denoted by ^a{x,t). In the time development of ^a, the control parameter e^ 
is not necessary to be constant. However in this article, we assume that e^ 
is sufficiently small for every t E T. 

Let the density of each ^^ be denoted by pa- We have the global density 
function p{x,t) and pressure pp{x,t) given by 

P(a;,i) = ^Pa^aix,t), pp{x,t) = ^Pa^aix,t). 

In contrast to the previous subsection, in this subsection, we investigate 
an initial problem. In other words, every configuration of the geometrical 
objects. Ma, La and approximately Ba (a = 0, ■ ■ ■ , A^ — 1), with divergence 
free velocity u, (div(M) = 0) can be an initial condition to the dynamics of 
the multi-phase fields. The following equations which we will derive in this 
subsection govern the deformations of these geometrical objects as their time- 
development. Further it is noticed that in this subsection, the proper pressure 
Pp{x, t) has no mathematical nor physical meaning because it becomes a part 
of the total pressure p, which is determined by the divergence free condition 
div(u) = as mentioned in Remark [H 

We have the first theorem; 

Theorem 1. The action integral of the multi-phase fields, or the e^-controlled 
color functions ^a with physical parameters pa, (Tab, Pa (a, 6 = 0, 1, ■ ■ ■ , A^ — 1) 
defined above, is given by 



5, 



mul 



I dt I (]-p\u\'' -Y,''-b^/\VU\Vib\{ia + ib) + pp] d^x, (5.11) 

'''^ "^^ V a>b J 
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under the volume-preserving deformation. 

Proof. The action integral is additive. The first term exhibits the kinematic 
energy of the fluids. The second term represents the surface energy up to 
e^ as in Proposition [121 The proper pressure pp in (15.51) leads the Laplace 
equations. We can regard it as the action integral of the multi-phase fields 
with these parameters. D D 



Then we have further generalization of (14.81) as follows: 

Lemma 2. Assume that every Ma{t), M^™^(t) and M^^°^(t) deform for the 
time-development following a certain equation. The Euler- Lagrange equation 
of the action integral with respect to C,a, S<Smui/S(,a = 0, is given, up to the 
volume preserving condition, as follows: 

1. For a point x G M^™^, we have 

idpa- Pb)\u{x,t)\'^ + (Pa-Pb) -aabKaix,t) = 0. (5.12) 

2. For a point x G M^^"^ , we have 

^{Pa - Pb - Pc)\u{x,t)\'^ + (Pa-Pb-Pc) - Kabc{x,t) = 0. (5.13) 

Similarly we have the similar equations for higher intersection regions. 

Proof. It is the same as proof of Proposition [131 D d 

Using these equations, we have the second theorem, which is our main 
theorem: 

Theorem 2. For every (x,t) G QxT, the variational principle, 6Saiui/Sl{x,t) 
0, provides the equation of motion. 



^diP + J2 ^a,b [di (V|Vea||V6|(ea + 6 



Dpu^ 

a>b 






(5.14) 



Here p is the pressure coming from the effect of the volume-preserving or 
incompressible condition, which includes the proper pressure pp ( 15.5]) . 
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Proof. We naturally obtain it by using 1) Proposition H] and its proof, 2) 
Remark [H 3) Lemma [2] and 4) Proposition El D D 

Here we note that by expressing the low-dimensional geometry in terms 
of the global smooth functions ^'s with finite e^, we have unified the infi- 
nite dimensional geometry or the incompressible fluid dynamics governed by 
IFluid(r2 X T), and the e^-parameterized low dimensional geometry with sin- 
gularities to obtain the extended Euler equation (15.141) . When e^ approaches 
to zero, we must consider the hyp erf unctions [261 129] instead of C°°{Q x T), 
but we conjecture that our results would be justified even under the limit; 
the unification would have more rigorous meanings. 

It should be noted that on the unification, it is very crucial that we express 
the low-dimensional geometry in terms of the global smooth functions ^'s as 
the infinite-dimensional vector spaces. The SDiff(r2) naturally acts on ,^'s 
and thus we could treat the low-dimensional geometry and the incompressible 
fluid dynamics in the framework of the infinite dimensional Lie group |H UHl 
[38] . It is contrast to the level-set method. As mentioned in Section [2Tn the 
level-set function does not belong to C°°{Q) and thus we can not consider 
SDiff (fi) action and treat it in the framework. 

Remark 11. 1. fl5.14p is the Euler equation with the surface tension to 
multi-phase fields which gives the equation of motion of the multi- phase 
flow even with the multiple junctions. As we will illustrate examples in 
Section El the dynamics with the triple junction can be solved without 
any geometrical constraints. It should also noted that for a point in 
M^j°^, (15.141) is reduced to the original Euler equation in Reference [31] 
or dHU). 

2. The Euler equation (I5.14p appears as the momentum conservation in 
the sense of Noether's theorem (Section 12. 2p . It implies that (I5.14p is 
natural from the geometrical viewpoint [H HI [TSl 1251 EBl ESI I 



3. Further even though we set {^a{-,t)} as proper e^-controUed colored 
functions as an initial state, their time-development is not guaranteed 
that {^a{-,t)}, (t > 0), is proper e^-controlled. In general e^ may be- 
come large for the time development, at least, numerically due to the 
numerical diffusion. (See examples in Section [6]). However even for 
t > 0, we can find e^(t) such that {C,a{-,t)} are eg(t)-controlled col- 
ored functions and if e^{t) is sufficiently small, our approximation is 
guaranteed by e^(t). 
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4. The surface tension is also defined over Q x T and thus the Euler 
equation is defined over QxT without any assumptions due to Remark 

m 

5. We may set e^ depending upon the individual intermediate region be- 
tween these fields by letting eat mean that for ^^ and ^t,, a ^ b. Then 

N-l 

if we recognize e^ as maxeat, above arguments are applicable for the 

a, 6=0 

case. 

6. We defined the e^-controlled colored functions using the er-tubular 
neighborhood Tjj^^r^ and the colored decomposition of fl in Definition H] 
by letting et = e^. On the other hand, as in Reference [31], our formu- 
lation can describe a topology change well following the Euler equation 
f l5.14p such as a split of a bubble into two bubbles in a liquid. The 
e^-controUed colored functions can represents the geometry for such 
a topology change without any difficulties. However on the topology 
change, the path-connected region and the e^-tubular neighborhood 
lose their mathematical meaning and thus, more rigorously, we should 
redefine the e^-controUed colored functions. Since the e^-controUed col- 
ored functions represent the geometry as an analytic geometry, it is 
not difficult to modify the definitions though it is too abstract. In 
other words, we should first define the e^-controlled colored functions 
^'s without the base geometry, and characterize geometrical objects us- 
ing the functions ^'s. However since such a way is too abstract to find 
these geometrical meanings, we avoided a needless confusion in these 
definitions and employed Definition HJ 

5.5 Equation of motion of triple-phase flow 

Let us concentrate ourselves on a triple-phase fiow problem, noting f l5.3p . 
From the symmetry of the triple phase, we introduce "proper" surface tension 
coefficients, 

C^Ol + O"02 - 0-12 O"01 + 0"12 - O'02 Cro2 + (T12 - <Toi 

^0 = ^ , a, = , a, = , 

or (Tab = cTa + cTb. Here it should be noted that the "proper" surface tension 
coefficient is based upon the speciality of the triple-phase and does not have 
more physical meaning than above definition. 
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Lemma 3. For different a, b, and c, we have the following approximation, 

(5.15) 
Using the relation, the free energy (15.11) has a simpler expression up to 

Proposition 15. By letting 

^sym := ^0 / d'x \Vio{x)\+ai / d^'x |Vei(x)| +a2 / t/^x |Ve2(x)|, 
Ju Jn Jn 

we have a certain number M related to area of the surfaces {Ba} such that 

Proof. Due to Lemma [3], it is obvious. D D 

The action integral (15. lip also becomes 

Stri = dt I 2^1"!^ ~ X^(cra|V^a| - PaL) ) d^X. 

For a practical reason, we consider a simpler expression by specifying the 
problem. 

5.6 Two-phase flow and wall with triple-junction 

More specially we consider the case that ^o corresponds to the wall which 

does not move. For the case, we can neglect the wall part of the equation, 

(s) 
because it causes a mere energy-shift of ^sym. Then the action integral and 

the Euler equation become simpler. We have the following theorem as a 

corollary. 

Theorem 3. The action integral of two-phase flow with wall is given by 

5wall = dt -P|MP - V'(cra|V^a| - Pa^a) d^X, 
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and the equation of motion is given by 

Dpu 



Dt 

where 



+ dip-dj{Ti,)=0, (5.16) 



^^"(^-^«^)i^^"i- ^''"^ 



a=l 



Practically this Euler equation fl5.16p is more convenient due to the proper 
surface tension coefficients. However this quite differs from the original (14. 5 p 
and (14. 6 p in Reference [31] and governs the motion of two-phase flow with a 
wall completely. 

Remark 12. Equation (I5.16P is the Euler equation with the surface ten- 
sion for two-phase fields with a wall or triple junctions in our theoretical 
framework. We should note that under the approximation ( JS.lSp . ( J5.16P is 
equivalent to (I5.14p . even though (I5.16P is far simpler than (I5.14p . 

From Remark [2], it should be noted that r and the Euler equation (I5.16P 
are defined over QxT. This property as a governing equation is very impor- 
tant for the computations to be stable, which is mentioned in Introduction. 
Since the non-trivial part of r is localized in Q of each t & T,t vanishes and 
has no effect on the equation in the other area. 

We will show some numerical computational results of this case in the 
following section. There we could also consider the viscous stress forces and 
the wall shear stress. 



6 Numerical computations 

In this section, we show some numerical computations of two-phase flow 
surrounded by a wall obeying the extended Euler equation in Theorem [3l 
As in Theorem [3l the wall is expressed by the color function ,^o and has the 
intermediate region (Mq \ Lq)° where ^o has its value (0,1). As dynamics 
of the incompressible two-phase flow with a static wall, we numerically solve 
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the equations, 



div(M) = 
Dpu' 



Dt 
Dp 

Dt ' 



0. 



Here for the numerical computations, we assume that the force K consists 
of the surface tension, the viscous stress forces, and the wall shear stress, 

Kj = difij + diTij + fj. (6.2) 

Here r is given by f l5.17p . r is the viscous tensor, 

/ 1 , A ^ 1 fdu' du^\ 

^- '■= 2^ [^- - 3^"("V ' ^- '■=2[d7, + d^J 

with the viscous constant 

r]{x) := 771^1 + 7/26, 

and fj is the wall shear stress which is localized at the intermediate region 
(Mo \ Lq)° where 6 has its value (0, 1). 

The boundary condition of the interface between the fluid 6 ('^ = 1, 2) 
and the wall C,o is generated dynamically in this case. In other words, in order 
that the wall shear stress term suppress the slip over the intermediate region 
{Mq\Lq)° asymptotically t — ?■ 00 due to damping, we let fj be proportional to 
j-component of du^^/dqo for the parallel velocity m" to the wall and relevant to 
{1 — C,o{x)) , and make u vanish over Lq. Here go, Mq, and Lq are of Definition 

il 

The viscous force can not be dealt with in the framework of the Hamil- 
tonian system because it has dissipation. However from the conventional 
consideration of the balance of the momentum [THl Sec. 13], it is not difficult 
to evaluate it. The viscosity basically makes the numerical computations 
stable. 

In the numerical computations, we consider the problem in the structure 
lattice C marked by aZ'^, where Z is the set of the integers and a is a positive 
number. The lattice consists of cells and faces of each cell. Let every cell be 
a cube with sides of the length a. We deal with a subspace Qc of the lattice 
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as f2£ := r2 n £ C E^. The fields ^'s are defined over the cells as cell wise 
constant functions and the velocity field u is defined over faces as facewise 
constant functions [15]; ^ is a constant function in each cell and depends on 
the position of the cell, and similarly the components of the velocity field, u^, 
M^, and u^ are facewise constant functions defined over x^x^-iaces, x^x^-iaces, 
and a:;^a;^-faces of each cell respectively. 

As we gave a comment in Remark [TT] 5. we make the parameter e^ depend 
on the intermediate region in this section. Let ei2 be the parameter for the 
two-phase field or the liquids, and eo := eoi = ^02 be one for the intermediate 
region (Mq \ -^0)° between hquids and the wall. 

As mentioned in Introduction, we assume that ei2 for the two-phase field 
in our method is given as ei2 > a so that we could estimate the intermediate 
effect in our model following References [3 HH HH |2H |3T], even though the 
thickness of the intermediate region among real liquids is of atomic order and 
is basically negligible in the macroscopic theory. 

In the computational fiuid dynamics, the VOF (volume of fiuid) method 
discovered by Hirt and his coauthors [151 Ell [22] is well-established when we 
deal with fiuid with a wall. Since we handled triple-junction problems as in 
Section 15. 6[ we reformulate our model in the VOF-method. It implies that 
we identify 1 — ^0 with the so-called V^- function V := 1 — ^0 in the VOF 
method because V in the VOF method means the volume fraction of the 
fiuid and corresponds to 1 — ^0 in our formulation. 

As the convention in Reference [22] , V is also defined as a cellwise constant 
function. In the following examples, we will set eo to be a or the unit cell 
basically. However we can also make it eo > a as for two-phase field. It 
means that for the case eo > a, we consider each cell as a fictitious porous 
material whose volume ratio V & [0,1] without imposing any wall shear stress 
on the fictitious surface of the porous parts itself in each cell as in Figure 
1. (As mentioned above, we set the wall shear stress fj from the physical 
wall ^Q. The porous parts are purely fictitious.) The region where V is equal 
to 1 means the region where fiuid freely exists whereas the region where V 
vanishes means the region where existence of fiuid is prohibited. The region 
with V G (0, 1) is the intermediate region (Mq \ Lq)°. Here we emphasize 
that the fictitious porous in each cell brings purely geometrical effects to this 
model. 

Then we could go on to consider the problem in consistency between 
VOF-method and ^0 function in the phase-field model. Let functions fi = f 
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and /2 over supp(K) be defined by the relations, 

^l = Vh, 6 = ^/2, /l + /2 = l. 

Further we also modify the open fraction A in the VOF-method, which 
is defined over each face. We interpret A as the open area of the fictitious 
porous material of each face of each cell, which also has a value in [0, 1] as 
in Figure 1. We also use the open area fraction A of each face of each cell 
[221 EI]- For a face belonging to the cell whose V = 1, A is also equal to 1. 
Following the convention in discretization by Hirt [22], A is regarded as an 
operator acting on the face-valued functions like 

Aou = Au = (Aim\ A2M^ AsU^), 

{AuY = Aiu\ {AuY = A2u\ {Auf = Asu^ (6.3) 

Here we note that AjO^ implicitly appearing in (16.31) can be interpreted as a 
two-chain of homological base associated with a face of a cell. For example, 
for a velocity field fi := u^{x)dx^ defined over a cell in the continuous theory 
and a piece of the boundary element of the cell Aio^, the discretized u^ 
defined over the face is given by 



{AuY := ^ / */i = Aiu\ 



where * is the Hodge star operator, i.e., *^ := u^{x)dx^dx^ + v?{x)dx^dx^ + 
v?{x)dx^dx'^. Thus the discretization (16. 3 p is very natural even from the 
point of view of the modern differential geometry. 

Hence div(ii) = Vu reads VAu as the difference equation in VOF-method 
[22] and we employ this discretization method. 

We give our algorithm to compute (16. ip precisely as follows. As a con- 
vention, we specify the quantities with "old" and "new" corresponding to 
the previous states and the next states at each time step respectively in 
the computation. In other words, we give the algorithm that we construct 
the next states using the previous data by regarding the current state as an 
intermediate state in the time step. We use the project-method [T3| [T5]: 

I : ^^^ = -(,.- . V)p.-. 

now _ ~ 1 

II:!4^^-i,Vp-/0. 
HI : Vm^^" = 0. 
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The step I is the part of the advection of the velocity u"^'^. In the step I, we 
define an intermediate velocity u and after then, we compute u'^'^'" and p in 
the steps II and III. 

The time-development of p is given by the equation, 

ynew _ yold ^ AtV((AM°'<^)/°''i), 

and 

P=V{p^f + p2{l-f)) 

for the proper densities pa of C,a (a = 1, 2). 

Even for the case that we can deal with multi-phase flow with large density 
difference, we evaluate its time-development. Precisely speaking, when we 
evaluate u, following the idea of Rudman [51] we employ the momentum 
advection u of u, 

1 



u ■= — — [p°i%°id _ At(M°^d . V)p°'%' 



Our derivation of the Euler equation shows that the Rudman's method is 
quite natural. 

Following the conventional notation, the guessed-value of the velocity is 
denoted by u*, which is the initial value for the steps in II and III. Let us 
define 

u* ■.= u + At^K(p°^'^, /°'^ u°^'^). 

In order to evaluate the guessed velocity, we compute the force K from (16.21) 
noting that divr and divr read VAt and VAt respectively. 

Following the SMAC (Simplified-Marker-and-Cell) method [31 [131 [15], 
we numerically determine the new velocity m"^*^" and the pressure j9 in a cer- 
tain boundary condition using the preconditioned conjugate gradient method 
(PCGM): 

(Ila) Evaluate p using the PCGM : -r-V(v4 ou*) = VAo ^^Vp, 

(lib) By using p determine u"™ : m'''^" = u* - At^^Vp. 

More precisely speaking, (III) V(v4 o m'^''") = means that we numerically 
solve the Poisson equation. 



[AAt— 



V AAt Vp = V(A o u* 
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Then we obtain m'^'^", which obviously satisfies (III) V{A o u^'^'^) = 0, which 
is known as the Hodge decomposition method [5], [121 [E] as mentioned in 
Remark [H 

Following the algorithm, we computed the two-phase flow with a wall and 
triple junctions. We illustrate two examples of the numerical solutions of the 
triple junction problems as follows. 



6.1 Example 1 

Here we show a computation of a capillary problem, or the meniscus oscil- 
lation, in Figure 2. We set two liquids in a parallel wall with the physical 
parameters; rji = rj2 = 0.1[cp], pi = p2 = l-0[pg//um^], cri = 3.349 [pg//isec^], 
(T2 = 46.651 [pg//isec^]. 

We used £ := 12[/im] x 0.5[/im] x 16[;um] lattice whose unit length a is 
0.125[/im]. The flrst liquid exists in the down side and the second liquid does 
in the upper side in the region 10[/im] x 0.5[/im] x 15[yum] surrounded by 
the wall and the boundaries with the boundary conditions. As the boundary 
conditions, at the upper side from the bottom of the wall by 15[yum], we flx 
the constant pressure as 100 [KPa] and, along a;^-direction, we set the periodic 
boundary condition. 

We set ei2 = eo = 1 mesh for the intermediate regions, at least, as its 
initial condition. Each time interval is 0.001 [/xsec]. 

As the initial state, we start the state that the fluid surface is flat as in 
Figure 2 (a) and the flrst liquid exists in the box region 10[/im] x 0.5[/im] x 
7.0[/im], which is not stable. Due to the surface tension, it moves and starts 
to oscillate but due to viscosity, the oscillation decays. Though we did not 
impose the contact angle as a geometrical constraint, the dynamics of the 
contact angle was calculated due to a balance between the kinematic energy 
and the potential energy or the surface energy. The oscillation converged to 
the stable shape with the proper contact angle, which is given by 

0-2 - O"! O"02 - O-Ql , - 

COS(y9 = = . (6.4) 

0-2 + CTl cri2 

The angle given by a's are designed as 30 [degree] whereas it in the numerical 
experiment in Figure 2 is a little bit larger than 30 [degree], though it is 
very difficult to determine it precisely. However since we could tune the 
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parameters cr's so that we obtain the required state, our formulation is very 
practical. 

Due to the numerical diffusions and others, the thickness of the inter- 
mediate regions changes in the time development and also depends on the 
positions of the interfaces, even though it is fixed the same at the initial 
state. However we consider that it is thin enough to evaluate the physical 
system since the contact angle is reasonably estimated. 

6.2 Example 2 

This example is on the computations of the contact angles for different surface 
tension coefficients displayed in Figure 3. 

Even in this case, in order to see the difference between the designed 
contact angle and computed one, we go on to handle two-dimensional sym- 
metrical problems though we used three-dimensional computational software. 
In other words, we set that x^-direction is periodic. 

Since the contact angle y? in our convention is given by the formula (16. 4p . 

By setting cr's 

(J I 1 — cos ip 

(72 1 + cos if ' 

for given the contact angle y?, we computed five triple junction problems 
without any geometrical constraints; each a is given in the caption in Figure 
3. The other physical parameters are given by i]i = ri2 = 0.1 [cp] and pi = 
p2 = 1.0[pg//im% 

In this computation we used a 240 x 4 x 112 lattice whose unit length a 
is 0.125[/im]; Q = 30[/im] x 0.5[/im] x 14[/im]. We set the fiat layer as a wall 
by thickness 3[/im] from the bottom of Q along the z-axis. As the boundary 
conditions, at the upper side from the bottom of the wall by 9[/im], we fix 
the constant pressure as 100 [KPa]. 

As the initial state for each computation, we set a semicylinder with 
radius 5[/im] in the fiat wall like Figure 3 (d). We also set ei2 = eo = 1 mesh 
for the intermediate regions. Each time step also corresponds to 0.001 [sec]. 

Due to the viscosity, after time passes sufficiently 50[/xsec], the static 
solutions were obtained as illustrated in Figure 3, which recover the contact 
angles under our approximation within good agreements. 
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7 Summary 

By exploring an incompressible fluid with a phase-fleld geometrically [H HI 
[T8| |25| |28| |32| |37] . we reformulated the expression of the surface tension 
for the two-phase flow found by Lafaurie, Nardone, Scardovelli, Zaleski and 
Zanetti [31] as a variational problem. We reproduced the Euler equation of 
two-phase flow (14.111) following the variational principle of the action integral 
(J42D in Proposition [JOl 

The new formulation along the line of the variational principle enabled 
us to extend (14. lip to that for the multi-phase (A^-phase, N > 2) flow. 
By extending (14. lip , we obtained the novel Euler equation (I5.14p with the 
surface tension of the multi-phase flelds in Theorem |2] from the action integral 
of Theorem [1] as the conservation of momentum in the sense of Noether's 
theorem. The variational principle for the inflnite dimensional system in the 
sense of References [H HJ [18] gives the equation of motion of multi-phase flow 
controlled by the small parameter e^ without any geometrical constraints and 
any difficulties for the singularities at multiple junctions. 

For the static case, we gave governing equations (15. 6p . (15. 7p and (I5.10p 
which generate the locally constant mean curvature surfaces with triple junc- 
tions by controlling a parameter e^ to avoid these singularities. As the solu- 
tions of (14. 4p has been studied well as the constant mean curvature surfaces 
for last two decades ^7\ HSl [201 HZ], our extended equations (15. 6p . (15. 7p and 
(I5.10p might shed new light on treatment of singularities of their extended 
surfaces, or a set of locally constant mean curvature surfaces. (Even though 
we need an interpretation of our scheme, for example, it can be applied to 
a soap film problem with triple junction.) It implies that our method might 
give a method of resolutions of singularities in the framework of analytic 
geometry. 

By specifying the problem of the multi-phase flow to the contact angle 
problems at triple junctions with a static wall, we obtained the simpler Euler 
equation (I5.16P in Theorem|3l Using the VOF method [221 |2T], we showed two 
examples of the numerical computations in Section El In our computational 
method, for given surface tension coefficients, the contact angle is automat- 
ically generated by the surface tension without any geometrical constraints 
and any difficulties for the singularities at triple junctions. The computa- 
tions were very stable. It means that the computations did not collapse nor 
behave wildly for every initial and the boundary conditions. 

In our theoretical framework, we have unified the infinite dimensional 
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geometry or an incompressible fluid dynamics governed by IFluid(r2 x T), 
and the e^-parameterized low dimensional geometry with singularities given 
by the multi-phase fields. We obtained all of equations following the same 
variational principle. We naturally reproduced the Laplace equations, (14. 4p 
and (15. 6p . and obtained their generalizations (14. Sp . (15. 6p . (15. 7p . (I5.13P and 
f lETU]) . and the Euler equations, fHTTT]) . f lETip . and flHTTB]) in Proposition HD] 
and Theorems [2] and O These equations are derived from the same action 
integrals by choosing the physical parameters. In the sense of References 
[H [5l [To], it implies that we gave geometrical interpretations of the multi- 
phase fiow. Even though the phase-field model has the artificial intermediate 
regions with unphysical thickness e^, our theory supplies a model which shows 
how to evaluate their effects on the surface tension forces, from geometrical 
viewpoints. The key fact of the model is that we express the low-dimensional 
geometry in terms of the infinite-dimensional vector spaces, or global func- 
tions ^ 's which have natural Diff and SDiff actions. Thus we can treat them 
in the framework of infinite dimensional Lie group |H UHl [38] to consider its 
Euler equation. It is contrast to the level-set method; in analytic geometry 
and algebraic geometry, zeros of a function expresses a geometrical object 
and thus the level-set method is so natural from the point of view. However 
as mentioned in Section 12. ![ the level-set function cannot be a global func- 
tions as C°°(i7) and thus it is difficult to handle the method in the framework 
of the infinite dimensional Lie group SDiff (^2). 

As our approach gives a resolution of the singularities by a parameter eg, 
in future we will explore topology changes, geometrical objects with singular- 
ities and so on, more concretely in our theoretical framework. When e^ ap- 
proaches to zero, we need more rigorous arguments in terms of hyperfunctions 
but we conjecture that our results would be correct for the vanishing limit 



of ec because the Heaviside function is expressed by 6{q) = lim — tan 

in the Sato hyperfunction theory, which could be basically identified with ^(g) 
of the finite e^. Since an application of the Sato hyperfunction theory to fiuid 
dynamics was reported by Imai on vortex layer and so on [29] , we believe that 
this approach might give another collaboration between pure mathematics 
and fiuid mechanics. 
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position 



Figure 1: VOF with porous matter expression: For the consistency between 
the color function method and VOF-method, we consider each cell as a fic- 
titious porous material whose volume ratio and open fraction are a value in 
[0, 1] without imposing any wall shear stress on fictitious surface of the porous 
parts in each cell. This expression represents purely geometrical effects. 
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Figure 2: The meniscus oscillation: Each figure shows the time development. 
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Figure 3: The different contact angles are illustrated due to the different 
surface energy: By fixing a\ = 1.0000[pg//isec^], (a): ip = 30 [degree], o"2 = 
13.9282[pg///sec3], (b): v? = 45 [degree], ^2 = 5.8284[pg//isec3], (c): ^ = 60 
[degree], (72 = 3.0000[pg//isec^], (d): ip = 90 [degree], (72 
and (e): ip = 120 [degree], (T2 = 0.3333[pg//isec^]. 
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